I =[];

I(1).pix = imread('band1.gif');
I(2).pix = imread('band2.gif');
I(3).pix = imread('band3.gif');
I(4).pix = imread('band4.gif');
[height,width]= size(I(1).pix);

for i= 1:height,
    for j=1:width,            
        data(1,i,j)= I(1).pix(i,j);            
        data(2,i,j)= I(2).pix(i,j);  
        data(3,i,j)= I(3).pix(i,j);  
        data(4,i,j)= I(4).pix(i,j);
    end % end for j
end % end for i
    
data=double(data);

dim= 4; 
N= height * width;
c= input( 'enter c (# of Classes) in Fuzzy_c_means algo :');
m= input( 'enter m ,the exponential weight :');
epsilon= input( 'enter epsilon, the termination criterion :');
flag= 1;
Normalize= [];
U0= rand(c,height,width);
Normalize= sum(U0,1);
for i= 1:height,
    for j= 1:width,
        U0(:,i,j)= U0(:,i,j)/ Normalize(1,i,j);
    end % end for j
end % end for i

U1= zeros(c,height,width);
v= [];
for i= 1 :c,
    v(i).mean= zeros(1,dim);
    v(i).denom= 1;
end % end for
% computation of new v(i)'s
while flag
 U0= U0.^m;
 for i= 1:c,
     S(i)= 0;
     for k= 1:height,
         for t= 1:width,
             S(i)= S(i)+ U0(i,k,t);
         end % end for t
     end % end for k
     v(i).denom= S(i);
     numerator= zeros(1,dim);
     for j= 1 : height,
         for k= 1: width,
             numerator= numerator + (U0(i,j,k)*data(:,j,k))';
         end % end for k
     end % end for j
     v(i).mean= numerator/(v(i).denom);
 end % end for i
 
 % computation for U1
 inv_norm= zeros(c,height,width);
 inv_pow= 1/(m-1);
 for i= 1:c,
     for j= 1: height,
         for k= 1: width,
             temp= (sum(data(:,j,k)- v(i).mean').^2)^inv_pow;
             if temp ~= 0,
                 inv_norm(i,j,k)= 1/temp;
             end %end if
         end % end for k
     end % end for j
 end % end for i
 
%while flag 
for j=1:height,
    for t= 1: width,
        denom= 0;
        for k= 1:c,
            denom= denom + inv_norm(k,j,t);
        end % end for k
        for i= 1:c, 
            if denom ~= 0,
                U1(i,j,t)= inv_norm(i,j,t)/denom ;
            end % end if
        end % end for i
    end % end for t
end %end for j

 if max(max(max( abs(U1- U0))))<= epsilon 
     flag= 0;
 end % end if
U0= U1;
end % end while

[Y,I]= max(U0,[],1);
image= zeros(height,width);
I= int16(I);
for i= 1:height,
    for j= 1:width,
        if I(1,i,j)== 1,
            image(i,j)=0;
        else if I(1,i,j)== 2,
                image(i,j)= uint8(70);
            else if I(1,i,j)== 3,
                    image(i,j)= uint8(140);
                else if I(1,i,j)== 4,
                        image(i,j)= uint8(210);
                    else image(i,j)== uint8(255);
                    end
                end
            end
        end
    end
end

    
%     C(j)=I(j);
% end
% plot(D,C,'o');



